package gates;

import master.*;
import java.math.*;

/**
 * 
 * @author Rupert
 * 
 * Shors factorization algorithm. To factorize an integer N this creates a random integer
 *  0 < a < N and computes the period of f(x) = a^x mod N by using a fourier transformation
 *  and contineous fractions. 
 *  
 *  The algorithm uses 2 QuantumRegister and simulates their entanglement. 
 *  The function f(x) = a^x mod N is implemented in a classical way, however functions of
 *  this type can efficiently computed on a quantum computer using basic gates.
 *  
 *  Shors fails for certain combinations of a and N. The algorithm here test if this is the
 *  case by making a GCD test at the end and reruns automatically if GCD = 1;
 * 
 *  If you try to factorize a prime it will run endless. 
 *
 *
 */

public class Shors extends QuantumGate {
	
	int N;
	int m;
	int a;
	int j;
	public Shors(int N) {
		this.N = N;
		this.m = (int)(Math.floor(Math.log(N*N)/Math.log(2))+2);
		this.j =(int)Math.pow(2, m);
	//	this.m = 13;
	}

	/** shors main routine, returns a fator of N
	 * 
	 * @return
	 */
	
	public int runShors() {
		
		a = (int) (Math.random() * ((N-1))+2);
		int result=0;
		int entangle[] = new int[j];
		QuantumRegister reg1 = new QuantumRegister(m);
		QuantumRegister reg2 = new QuantumRegister(m);
		
		// create a equal superposition of all states in reg 1
		reg1.setToEigenstate(0);
		QuantumGate gate = new FullHadamard();
		gate.actOnReg(reg1);
		
		
	/* This is a shortcut routine that prevents the actual algorithm from beein executed when a 
		already is a factor of N;
		if(GCD(a,N)!= 1 && (a!=N)) { return GCD(a,N); } 
	*/	
		entangle = functionX(reg1,reg2);
		
//		measure 2nd register (f(x)).
		reg2.collapseToEigenstate();
		
	/*	find value in 2nd register, this could be more efficient if collapseToEigenstate() would return
		the eigenstate it collapsed to. This is however closer to the actual QM behavior 
	*/
		for(int n=0;n<j;n++) { if(reg2.getCoefficient(n).getReal()==1) result = n; }
	
	/*  simulated the entanglement of the first register to the measurement in the second and renormalize.
	 */
		for(int n=0;n<j;n++) { if(result != entangle[n]) reg1.setCoefficient(new ComplexNumber(0,0), n); }
		reg1.normalise();
		
//		apply fourier to reg1 and measure again
		QuantumFourier.applyFourier(reg1);
		reg1.collapseToEigenstate();
		for(int n=0;n<j;n++) { if(reg1.getCoefficient(n).getReal()!=0) result = n; }
		
// 		find period and factor, test and return/run again		
		float z = (float)result/N;
		int r;
		
		r = denominator(z,(int)Math.pow(2, m/2));
		if(r%2!=0) { r = GCD(N,(int)(Math.pow(a,r/2)+1));}
		if(r == 1 || r == N) {return this.runShors(); }
		System.out.println("Shors succes");
		
		return r;
	}
	
	/**
	 * function that computes a^x mod N in an efficient way.
	 * @param a
	 * @param x
	 * @param N
	 * @return
	 */
	
	
	public static int functionmodN(int a, int x, int N) {
		
		if(x==0) return 1;
		int ret=a%N;
		int l = (int)Math.floor(Math.log(x)/Math.log(2));
		for (int n = 1; n<=l;n++){
			ret = (ret*ret)%N;
		}
		for(int n = 1; n<=(int)(x-Math.pow(2,l));n++) {
			ret = (ret *a)%N;
		}
		return ret;
	}
	
	/**
	 * 	 * This function calculates f(x) = a^x mod N for the values of the first register @param reg1
	 * and puts the result in the second register @param reg2. It returns an integer vector witch contains
	 * the entanglement information for the 2 Quantum registers.
	 * This implementation is purely classical, but simulates how an actual Quantum computer would
	 * work.
	 * @param reg1
	 * @param reg2
	 * @return
	 */
	
	public int[] functionX(QuantumRegister reg1, QuantumRegister reg2) {
		
		int fx;
		int entangle[]= new int[j];
		float amplitude[] = new float[reg1.getDimension()];
		for(int n = 0;n<j;n++) { amplitude[n]=0;}
		
		for(int n = 0;n<j;n++) {
			fx = functionmodN(a,n,N);
			amplitude[fx] = amplitude[fx] + reg1.magSquared(n).magnitudeFloat();
			entangle[n]=fx;
		}
		for(int n = 0;n<j;n++) {
			reg2.setCoefficient(new ComplexNumber((float)Math.sqrt(amplitude[n]),0f),n);
		}
		return entangle;
	}
	
	/** 
	 * Computes the greates common divisor using the GCD method of Math.BigInteger.
	 * @param a
	 * @param N
	 * @return
	 */
	
	public int GCD(int a, int N) {
		   BigInteger b1 = new BigInteger(""+a); 
		   BigInteger b2 = new BigInteger(""+N);
		   BigInteger gcd = b1.gcd(b2);
		   return gcd.intValue();
	}
	
	/**
	 * finds the denominator of the best rational approximation q2/q1 for x with the constraint
	 * q1,q2 < qmax. Uses contineous fraction expansion.
	 * @param x
	 * @param qmax
	 * @return
	 */
	
	public int denominator(float x,int qmax) {
		  float y=x;
		  float z;
		  int q0=0;
		  int q1=1;
		  int q2;
		  
		  while(true) {
		    z=y-(float)Math.floor(y);
		    if(z<0.5/(qmax*qmax)) { return q1; }
		    y=1/z;
		    q2=(int)Math.floor(y)*q1+q0;
		    if (q2>=qmax) { return q1; }
		    q0=q1; q1=q2;
		  }
	}
}
